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Abstract 

We study QCD at non-zero quark density, zero temperature, infinite cou- 
pling using the Glasgow algorithm. An improved complex zero analysis gives 
a critical point \x c in agreement with that of chiral symmetry restoration 
computed with strong coupling expansions, and monomer-dimer simulations. 
We observe, however, two unphysical critical points: the onset for the num- 
ber density jx , and fj, s the saturation threshold, coincident with pathological 
onsets observed in past quenched QCD calculations. An analysis of the prob- 
ability distributions for particle number supports our physical interpretation 
of the critical point [i c , and offers a new intepretation of fi Q , which confirms 
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its unphysical nature. The perspectives for future lattice QCD calculations 
of the properties of dense baryonic matter are briefly discussed. 
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I. INTRODUCTION 



Numerical simulations of lattice QCD at high temperatures are making quantitative 
theoretical predictions which will be confronted with experiments at RHIC and LHC [|l|. 
However, numerical simulations of QCD in an environment rich in baryons lags far behind. 
Phenomenologically we know that nuclear matter can exist up to a density of four times or- 
dinary matter in neutron stars, and that higher density will eventually induce deconfmement 
and chiral symmetry restoration because of asymptoptic freedom. Current estimates from 
phenomenological nuclear models || place the critical chemical potential between 1000 and 
1600 Mev, and the critical baryon density between twice or twenty times that of ordinary 
nuclear matter. 

As it is well known, the reason behind this poor knowledge is the lack of a reliable calcu- 
lational scheme for lattice QCD at high baryon densities 0. A solid theoretical formulation 
for finite density QCD was made ten years ago ||, 0, but, since the resulting action is 
complex, probabilistic simulation methods fail. Early approaches considered the quenched 
approximation, which omits the complex fermion determinant, but it produced unphysical 



results [pi- [101- To obtain reliable results the determinant should be included, and the 



Glasgow method has been proposed to tackle this challenge ||11|| - fllTf . Although we are 
ultimately interested in the weak coupling, continuum limit, the strong coupling limit of 
QCD is attractive for several reasons : 1. There are analytic results coming from the strong 
coupling expansion []18 |, [19| , |p0l and numerical results from monomer-dimer simulations 
PHI , 2. The theory confines and spontaneously breaks chiral symmetry. In this paper we 
will use these features of strongly coupled lattice QCD to test and shed light on simulation 
methods which could be used at any coupling |22| . 

This paper is organized into two main Sections, Method and Results. 

Method is part review, and part illustration of our method of analysis. We first review 
the Glasgow method, and the relevant observables (II. A). We continue by illustrating some 
features which will help our numerical analysis. We will first discuss (II. B) the pathologies 
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found in calculations on isolated configurations. In Section II. C we will discuss how the 
Glasgow algorithm can escape from these single configuration pathologies, and build the 
physical signatures for the critical point /i c where chiral symmetry is restored. 

In Results, after discussing some generalities of the generation of the gauge ensemble 
(III. A), we show that the Glasgow method results for the Baryon current inherit some 
of the quenched/single configuration pathologies (III.B). Nevertheless, we will succesfully 
measure the critical chemical potential ji c (III.C) and we will discuss in detail the interplay 
of successes and failures of the results. Comprehensive results for an extended range of 
masses are given in Section III.D. In Section III.E we will present an alternative reanalysis 
of the critical region which will use probability distributions. We will confirm there the 
estimate of the critical point /i c , and we will offer a new intepretation of the pathological 
region. 

We conclude with a brief summary and discussion. 



II. METHOD 

Given the failure of the quenched approximation to deal with the problem of the chiral 

phase transition at high quark density, the natural conclusion is that dynamical quark 

simulations are essential. However, the complex measure of the functional integration with 

nonzero chemical potential poses a severe problem for such simulations. One simulation 

method which circumvents this problem is based on the expansion of the Grand- Canonical 

Partition Function (GCPF) in powers of the fugacity. The GCPF, (Z), can be written as 

the ensemble average ( |M(o'm)| ) wnere \M(fj,, m)\ is the fermion determinant at chemical 

potential \i and quark mass m in lattice units (lattice spacing a = 1), i.e. 

^ _ J[dU][dW}\M(^m)\e~ s ^ 

J[dU][dW}\M(0,m)\e- s ^ u ^ 1 } 

where U are matrices representing the gauge degrees of freedom and S g is the standard 
Wilson gauge action. The fermion matrix M is that describing 4-flavours of staggered 
fermions. 
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The fermionic determinant can be expressed explicitly as a function of \i by 

det(M(/x, m)) = e - 3 ^ nt det(P - e M ) (2) 
The lattice size is n^n t and P is the propagator matrix (independent of fi) [[|J 



' —GV V ^ 
y-V Oj 



(3) 



where G contains all the spacelike gauge links and the quark bare mass, and V all the 
forward timelike links of the fermionic matrix M. 

det(M(/x, m) can be computed in a basis where the propagator matrix is diagonal 

6nf rat 

det(M(//, m)) = e" 3 ^ nt n ( A fc ~ e ") ( 4 ) 

fc=i 

We recognize that the zeros of the determinant in the e M plane are the eigenvalues of the 
propagator matrix. The symmetry of the eigenvalues of the propagator matrix \k+j = 
e i2nj/n t ^ k f or j = o to n t — 1, together with the polynomial decomposition 

rat— 1 

n {e i2nj/nt (3 -x) = {f3 nt - x nt ) (5) 

3=0 

yields the equivalent representation 

6ra3 

det (M(p, to)) = e" 3 ^ ni ]J (A£* - eT*) (6) 

fc=i 

and dictates the general structure of the characteristic polynomial det(M(/j, to)) 

3raf 

det(M(/i,m))= ]T h i^ nt (7) 

k=— 3n| 

Note the dependence on /x is now via the fugacity / = e^ nt . 

Hence, measurement of the average of the characteristic polynomials (normalised by 
|M(0, to) |) in the ensemble generated at update mass to and ji = will give Z(fi, m) explicitly 
as a function of \x at that mass. 

This representation leads to a polynomial expansion of Z(fi) in powers of the fugacity 
whose coefficients are functions of the gluonic fields. 



3n 3 a 3n 3 3 

Z(n)= £ <k>e^ T = £ Q k f k (8) 

k=— 3n 3 k=— 3n| 

This expansion is just that of the GCPF expanded in terms of the canonical partition 
functions, (CPF's), for a fixed number of quarks (anti-quarks) on the lattice. Thermody- 
namical averages, which can be calculated as logarithmic derivatives of the GCPF, are then 
given explicitly as functions of fi. 

The relative value of the CPFs can characterise the properties of the system as well. 
For example, the relative weight of the triality-bearing to the triality-zero CPFs can signal 
whether the system is in the confined or deconfined phase. In the confined phase the ensemble 
average of the triality-bearing CPFs must be zero. This leads to 

n 3 

= E Qskf k (9) 

fc=n 3 

One can also explore the phase structure of the simulated system by examining the 
distribution of the zeros of the GCPF in the complex chemical potential (or fugacity) plane 
| 23fl 1 24]. These zeros correspond to the singularities of the thermodynamic potential and 
will converge in the thermodynamic limit, (L — > oo) towards any critical fi in the physical 
domain. 

In the following we also show that one can regard the zeros of the averaged characteristic 
polynomial as the "proper" ensemble average of the eigenvalues of the propagator matrix. 
This interpretation of the zeros as reflecting the ensemble average of the eigenvalues could 
be important in the interpretation of "unexpected" onset chemical potentials in ensembles 
of limited statistics. 



A. Observables 

The starting point of our analysis is the GCPF Z computed with the Glasgow algorithm. 
Our raw data are the CPF's Qn and our basic observables are the particle number density 
and the zeros of the GCPF in the complex fugacity plane. 
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Most of our discussions will consider the current < Jo(/i, to) > or equivalently the particle 
number density, defined as 

1 01n(Z(//,m)) _ 1 dln<det(MQu,m)) > 
< J (/i,m) >- - - - - - (10) 

Singular behaviour of the current can result from singularities in the density of baryonic 
states (particularly apparent in the zero-temperature limit). These singularities could be 
purely lattice artefacts and vanish in the continuum limit. However, they may instead 
reflect continuum spectral features, such as gaps in the spectrum or abrupt changes in 
the dispersion relation of the baryonic excitations. A chiral phase transition is one such 
possibility. A spectrum of chirally symmetric baryonic excitations will follow a gap less 
relativistic dispersion relation, contrary to the dispersion of particles with broken chiral 
symmetry. If the disappearance of the mass gap occurs together with the deconfinement 
transition, quark states will emerge instead of collective colourless baryonic excitations. 
Thus, the \i— dependence of J should determine the phase structure of dense baryonic 
matter, an alternative to the evaluation of the chiral condensate. 

Differentiating the action with respect to /i reveals the operator form of the charge, and 
one sees that the current is the expectation value of the number of paths through the links 
in the time direction 0] . In this sense the current can be defined on isolated configurations, 
where it reduces to 

4M = ^ dHiet( ^' mm (id 

In the quenched ensemble ln(det(M)) is differentiated before taking the statistical average: 

<, 1>Ifcm) = I/^!M\ (12) 

V \ dfi I 

and we recognize that 

< J (//, m) > q =< Jq(/jl, m) > (13) 

In the following the /i and m dependence will be left implicit wherever this does not create 
ambiguities. 



B. Failures on isolated configurations, and the quenched model 



The early work by Gibbs made it clear that the behaviour of some observables mea- 
sured on isolated configurations at finite density can be pathological. Since the analysis of 
isolated configurations is a necessary step in any lattice simulation, the impact of his result 
may be broader that its original motivation- to understand the pathologies of the quenched 
approximation. 

Our renewed interest was prompted by two considerations. First, our results presented 
in Section III. A below show clear relics of the quenched pathologies discussed in Gibbs's 
paper: the onset \x Q where the current Jo departs from zero is at half the pion mass. Second, 
published results on four fermion models |26| do not have such pathologies. However, 
both models share the same pattern of chiral symmetry breaking, and both models have 
Goldstone modes. Why, then, is there a difference at finite density? We decided to re- 
examine the behaviour of observables on isolated configurations in order to test the Gibbs 
scenario in a more general framework, and to gain some understanding of the process of 
statistical averaging in the two models. This paper is devoted to QCD, and the results for 
four fermion models will be presented elsewhere. 

First consider the behaviour of the current on isolated configurations. Jq follows from 
eqs. 0, | 

■>0 = -l + V (14) 

(here and in the following we use Gibbs' notation z = e M @). In the zero temperature case 
the sum over complex poles can be conveniently done by contour integration, yielding: 

Jl = ^7 E 1- (15) 

K|A 4 |<eM 

The threshold for the current Jq on isolated configurations is triggered by the lowest zero 
of the determinant. In turn, the zeros of the determinant are given by the eigenvalues of 
the Propagator Matrix that, as emphasized by Gibbs, are controlled by the mass spectrum 
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of the theory. The argument, which we briefly summarize for the sake of completeness, 
requires the calculation of the hadronic spectrum on replicated lattices, i.e. lattices which 
have been strung together d times in the time direction and the limit d — > oo is taken, in 
order to replace finite sums with contour integrations. This procedure is justifiable at zero 
temperature. 

The expression for the inverse of the fermion matrix, G(ti,t2), on the replicated lattices 
reads (slightly simplifying Gibbs' notation), 

G{t l M) = T. A ^ t r t2 (is) 

k 

where the A a are the amplitudes which can be related to the eigenvectors of the propagator 
matrix, and the are the corresponding eigenvalues. 

Eq. [16] shows that the exponential decay of G(ti,t 2 ) at large time-like separation is 
controlled by the eigenvalues of the propagator matrix. In other words, the eigenvalue 
spectrum calculated on isolated configurations should be closely related to the physical 
mass spectrum. In particular, Gibbs concluded that the smallest mass state is related 
to the lowest eigenvalue : 

= 21n|A min | (17) 

This identification was clear in simulations done by Gibbs because the pion propagator was 
very similar configuration by configuration although, strictly speaking, masses are properties 
only of the statistical ensemble. 

Gibbs argument has been reformulated and verified by Davies and Klepfish || . Patholo- 
gies of isolated configurations, the role of confinement and other issues, are also discussed 
in 0, [[K]]. All of these works confirm that on isolated configurations there is a singularity 
at a value of the chemical potential close to half the pion mass. 

Therefore, the results on isolated configurations are qualitatively different from those 
expected of the statistical ensemble. 

Some of the problems with the quenched model can be understood from eq. [H| : the 
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quenched current is a simple average of the one-configuration-current, and the quenched 
ensemble retains the pathological features observed on isolated configurations. 



C. The statistical ensemble, and the full model 

We can now focus on the interplay between the one configuration/quenched results and 
ensemble results. How can statistical averaging remove the problems observed on isolated 
configurations? Equivalently, how can the Glasgow averaging discussed above improve upon 
the quenched approximation? 

Consider the fugacity expansion for Z, eq. || By reinstating a factor e 3n ° nt ^ we see that 
Z is a polynomial of degree &n 3 n t = 6V in the variable z = e M . Z can then be written in 
terms of its zeros a, in the z plane 

Z = e 3V »Y[(z- ai ). (18) 

i=l 

Recall that Z =< det(M) > and compare formulae [18| and |] which we rewrite here 

6V 

detM = e 3 ^n( 2 -^) ( 19 ) 
i=i 

We see that the zeros of the partition function are the "proper" ensemble average of the 
eigenvalues of the fermionic propagator matrix, or, equivalently, of the zeros of the determi- 
nant. 

Manipulations analogous to those of eqs. 113, [T3| lead to the current 



■Jo = ^ E 1- (20) 

K|osj|<eM 



Let's search for other critical points past the first onset. From eq. ^0] we see that 
discontinuities in J are associated with a high density of zeros on circles with radius e Me in 
the fugacity plane. More generally, the density of the modulos of zeros in the e M plane is the 
derivative of Jq with respect to /i, i.e. the quark number susceptibility. Interestingly, the 
relevant quantities controlling the critical behaviour of the current are indeed the modulos 
of the complex zeros. 
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It is worth noticing that once the Z3 symmetry is enforced (eq. |^) 

2V 

Z = e 3n ^Y[(z 3 -p i ). (21) 
z=i 

The zeros in the complex plane z should then come in triplets, corresponding to cubic roots 
of certain complex numbers In principle (in practice things can be very different!) the 
effect of the Z3 symmetry can simply amount to a redistribution of phases with no effect 
on the moduli. That would not affect the critical behaviour, since the critical behaviour is 
triggered by the moduli themselves. The unphysical quenched onsets could certainly survive 
the Z3 symmetry of the full ensemble. 

The zeros of the partition function drive the critical behaviour of the full model as the 
zeros of the determinant drive the critical behaviour of isolated configurations, hence of 
the quenched model. In the process of going from the zeros of the determinant to the 
zeros of the grand canonical partition function, the pathological results observed on isolated 
configurations should turn into the physics of the full model : the fake critical points should 
disappear, the real phase transitions of the full model should emerge. 

III. RESULTS 

We present here our numerical results. All the background material, when not explicitly 
referenced, can be found in the previous Section. 

The generation of the configurations of gauge fields is described in The ensemble 
In The number density we review past results from the quenched approximation, from 
the monomer-dimer simulation of the full four flavor model and from the strong coupling 
expansion of the four flavor model. We present the results obtained with the Glasgow method 
on various lattice sizes, and we highlight the main similiarities and differences among the 
Glasgow results and the above-mentioned ones. A discussion of finite size effects is presented 
as well. In this and in the subsequent subsection, the emphasis is on the presentation of the 
main features of the results. We then limit ourselves to a discussion of one representative 
mass value, m q = .1. 

9 



In The determination of the critical point we focus on the analysis of the complex zeros 
in the e M plane. We contrast the pattern of zeros with that of the eigenvalues of the fermion 
propagator. We show how simulations of the full model produce a clear signature for the 
critical point fi c . 

In Light and heavy masses we will present our complete set of results. The light masses 
will offer information about the chiral limit. We confirm that our estimate // c remains 
constant, and different from zero when m — > 0. We demonstrate the dependence of \x c on 
the bare mass in the heavy quark regime. In the summary plot the results are compared 
with the strong coupling expansion, and the monomer-dimer calculations. 

The analysis of the probability distribution consists of a self-consistent analysis of the 
critical region which exploits the form of the probability distribution as a function of chem- 
ical potential. This analysis further validates our estimate of the critical point /i c and, in 
addition, offers a new intepretation of the onset region fi Q . 



A. The ensemble 

The Glasgow algorithm takes as input an ensemble of configurations at zero chemical 
potential. At infinite gauge coupling we can generate configurations either with the usual 
hybrid MonteCarlo procedure or just choosing random SU(3) matrices - this corresponds 
to different normalization for the partition function. First, we reproduced the results of 
Barbour, Davies and Sabeur |l]J , which show that the reweighting actually works on a 2 4 
lattice provided the statistics are high enough. Some preliminary runs were performed on 
a 4 4 lattice to check that the results were indeed independent of the algorithm chosen for 
the generation of the configuration. We finally selected a random generator which produces 
decorrelated configurations. 

We will present results on a 6 4 lattice for bare mass values ranging from .05 to 1.5 and 
on a 8 4 lattice for masses .08 and .1. The number of gauge field configurations analyzed 
ranges from a small sample of 25 on the 8 4 lattice m q = .08, ~ 100 configurations on the 
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same lattice, m q = .1, and several hundred configurations on the 6 4 lattices. 



B. The number density 

We studied m q = .1 on 6 4 lattices where we can contrast the results with those obtained 
1. in the quenched case, 2. with the monomer-dimer simulation, and 3. with the analytic 
results of the strong coupling expansions. Let's briefly review methods 2. and 3. 

The monomer-dimer approach (valid only at infinite coupling) writes the strong cou- 
pling action in a fashion suitable for computer simulations. It begins with the standard 
lattice QCD action with four flavors of staggered fermions and integrates out the completely 
disordered gauge fields. Confinement is enforced exactly and the short-ranged interactions 
between fermions allow the Grassman integrals to also be done exactly. The resulting action 
can be interpreted graphically in terms of 'monomers' and 'dimers', familiar constructions in 
statistical physics. This representation of the theory is well suited for computer simulations 
since the dreaded sign problem of the fermion determinant is not numerically significant in 
this basis (on small lattices). Computer simulations at m q = .1 [|2TJ] show a sharp transition 
at fi = .69(1) where the chiral condensate falls from its zero-/! value (essentially) to zero, 
and the number density J jumps from zero (essentially) to a fully occupied lattice, J = 1.0. 
These results agree with those of the traditional strong coupling expansion, as they should. 



The strong coupling expansion |T8[ |T9[ |20] at m q = .1 predicts, in fact, a strong first 



order transition at fi = .65 ( The small difference in fi c , can probably be accounted for by 
1/d corrections ). The analytic expressions of the strong coupling expansion show a feature 
not seen in the monomer-dimer simulations: a mixed phase for [i Q < < /i s where ordinary 
confined hadronic matter coexists with the saturated lattice phase. 

The quenched results [|Hj were characterized by a "forbidden region" ranging from fi Q = 
= . 32 to fi s ~m#/3~ 1.0. \i a and /i s are close to the extrema of the mixed phase 
predicted by the strong coupling expansions mentioned above. There is no remnant of the 
critical point for chiral symmetry restoration fi c ~ .65 predicted by the same expansion. 
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Our motivation in undertaking the 6 4 calculations, was, of course, to see results com- 
pletely different from the quenched calculations and very similar to the monomer-dimer 
results. 

These expectations were only partially borne out. We indeed found a signature at /x c , 
but also the persistence of \l q and \l s . These results are shown in Fig. [I], for the number 
density obtained with the Glasgow method on a 6 4 lattice at m q = .1. The Glasgow results 
are distinguished from the quenched ones by a small jump at \i ~ .7 ~ /i c , suggesting 
restoration of chiral symmetry. Other than that, the Glasgow and quenched results are very 
similar. 

We then moved to a larger lattice to study the sensitivity to size and temperature. Would 
the small hint of a discontinuity at fi ~ .7 become more pronounced? Would the dynamical 
results differ more substantially from the quenched ones? 

The answer was in the negative. Neverthless, we did learn something from these runs. 

In Fig. |2|, we show a detailed comparison of the results on the two lattices, 6 4 and 8 4 , 
m q = .1 (note the different scales on the right and the left side). By blowing up the picture 
of the number density, we see that the density itself deviates from zero at fi ~ 0.. This effect 
is very small (note the scale) and it is sensitive to temperature, the number density being 
suppressed, as expected, on the colder lattice. 

The most interesting point is that temperature effects are greatly lessened for fi > fi . 
One might have well thought that the increase at /i Q reflects a thermal excitation of baryons. 
This does not seem to be the case: only for /i < fi Q do we observe the expected, physical 
pattern of finite temperature effects. This disappears for /i > jj, . This result supports the 
belief that the rise at /i G is unphysical, as in the quenched approximation. Of course we 
cannot rule out the possibility that the situation changes on larger lattices, and we refer to 
[K| for discussions on this point. 

Temperature effects become apparent again at fi c , suggesting that fi c is a threshold of a 
new phase. 
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C. The determination of the critical point 



We can substantiate this intepretation of /x c by examining the zero's of the grand par- 
tition function in the complex plane e M . 

The numerical strategy suggested by Section II, eq. |2(] is straightforward: observe the 
distribution of the modulos of the zeros, or, equivalently, search for a strip of high densities 
in the e M plane. This criterion is numerically more convenient than the conventional Lee- 
Yang analysis, which only uses the zero whose imaginary part is closest to the real axes. 
It is also very natural : it says that the number density counts the density of states in the 
fermionic sector. 

This strategy is demonstrated in Fig. [| where we show the distribution of zeros (in 
practice, of the logarithms of their moduli) accompanying Fig. |l| The signal at \l c = .687(15) 
is very clear, and in excellent agreement with the monomer-dimer results fi c = .69(1). 

Figs. |] contrasts zeros of the determinant and zeros of the GCPF on a 8 4 lattice. The 
signal at /z c in the full model is quite clear. 

As discussed above, the upper and lower parts of the figure can also be seen as full 
and quenched results. As a by-product of our investigation, we see clearly why the search 
of further critical points in the quenched calculations was futile [^7| flTUfl : the eigenvalue 
distribution is almost "flat" . 

Unfortunately, the upper and lower parts of the figure also show strong signals at /i and 
ji s : the behaviour at the two side peaks does not change as we pass from the quenched to the 
full model. This gives us another view of the puzzling persistence of the onset noticed in the 
previous subsection. The Glasgow simulation method has failed to reproduce the published 
monomer-dimer results. 

We can also look at the zeros themselves, which we display in Fig. [5[ As anticipated in 
the discussion, we observe a dense line, which follows the prediction |e M | = e Mc . We also see 
a ring of zeros at half the pion mass, and note that the zeros fill up the entire region /i -yU s . 

In summary, we have seen how the small discontinuity observed in the current manifests 
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itself in the histogram of zeros: the density of zeros is the derivative of the number density, 
so a small "discontinuity" in J corresponds to a distinct signal in the histogram of zeros. 

Will more statistics eventually cancel the onsets at fi a and Even if we haven't ob- 
served any dramatic effect by increasing the number of configurations, insufficient statistics 
remains a possibility, especially since Z 3 invariance has not been completely achieved yet, 
and since it is possible that the precision required to achieve the cancellation of the un- 
wanted onset is prohibitively high. It is also possible that the polynomial representation for 
the GCPF is ill-conditioned |28[. 

D. Heavy and light masses 

All the results we have discussed above were for m q — .1. It is instructive to explore both 
light and heavy masses. Light masses are important for the chiral limit. There we expect the 
critical point fi c to remain constant and different from zero. Heavy masses change the critical 
point and allow a more detailed comparison of our results with the monomer-dimer /strong 
coupling approaches. 

Figs. ^ , [7| show the sensitivity to the quark mass close to the chiral limit : note the 
stability of the central peak, and the shift of the lower peak, which follows the onset of the 
current (see again also Fig. [|). 

Finally, we have also simulated larger masses, and the results are displayed in Fig. || We 
can appreciate the shifting in the central peak, and also its broadening - probably due to the 
fact that at large quark mass the transition is washed out. Interestingly, the current onset 
at larger quark mass is, apparently, smaller than half the pion mass - a surprising result 
since at m q = 1.5 in the quenched model the critical region shrinks to zero || - certainly 
this adds to the complication and the confusion associated with /i D . 

However, even if the interpretation of the critical region at this stage is largely subjective, 
the estimate of the critical point [i c seems reasonably sound. 

We then conclude our results Section with the summary of Fig. |9|. 



14 



E. The analysis of the probability distributions 



In this Section we re-examine the critical region by studying the probability distribution 
for the particle number. 
Write 

3V 

Z= £ W n (22) 

n=-3V 

and normalize such that Z = 1. W n is then the probability that a system in a grand 
canonical ensemble has n particles. 

Using the numerical results for the GPF above (see eq.(|8|), the shapes of the probability 
distributions W n = Q n e fintn for different chemical potentials can be drawn as a function of 
n, and the critical region can be studied using standard statistical mechanics analysis. 

For a transition in a classical ensemble in the infinite volume limit, the distribution of 
W n should have a single peak in a pure phase, and a flat distribution at the critical point, 
where all values of the particle number between the two extrema should be equally likely. 

The current < J > can be written as 

^ 3V 

<J o>=T7 E nW n (23) 

V n=-3V 



We draw in Fig. |10| the probability distribution for small chemical potential. At /i = 0. 
(solid line) the distribution is symmetric around the origin : look at the two satellite peaks 
of equal height. < J > equals zero as it should. At \i = .1 (dashed line) the distribution 
becomes asymmetric, reflecting the enhancement (suppression) of the forward (backward) 
propagation : the peak on the left decreases, the one on the right increases. Positive and 
negative states are still both contributing to the probability distribution. The net < J > 
moves immediately off zero, but it is very, very small (look again at the picture ^ ). The 
distribution broadens on smaller lattices, which accounts for the pattern of finite size effects 
seen in the same picture. 

At [L — /i the scenario changes completely: a secondary maximum develops at positive 
n, and the distribution moves to the positive n region. We show this behaviour for both 
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m q = .1 and m q = .08 in Fig. |TT[ For // > \i the negative states do not contribute. This 
behavior is correlated with the sharp increase of Jo plotted before and should be related to 
changes in the theory's spectrum, perhaps reflecting pathologies of the quenched case such 
as the " funny pions" |I0] or Stephanov's condensates ||29|| . The distribution is now roughly 
symmetric, and its broadening on smaller lattices does not affect its average value < J >. 
This behavior is compatible with the absence of strong volume effects, Fig.||]. 

Next, the critical region : we see the expected broadening of the probability distribution 
at (i c (Fig. , pD. Finite size effects become important again for fi > /i c . Note, in 



particular, in Fig. [13], the sensitivity to ji on a very fine scale : the three central plots are 
for fi = .68, .683, .7. 



In Fig. [14] we summarize these observations by plotting Wq, the probability that the 
system has zero particle number, and the integrated probabilities W + = J2W n ,n > 0; 
W~ = J2W n ,n < 0. The logarithmic scale of the plot makes it easy to see that backward 
and forward propagations are enhanced and suppressed by the same factor at small fi. 
Correspondingly, the contribution of n = must decrease. At // = [l n = equals the 
overall contribution from > 0. For \x > fi Q only positive n contribute to Z. 

These results suggest that \i Q is the threshold for a phase with only positive propagation. 
Perhaps this observation is a clue to the nature of the phase /i > \l q . Recall that mean 
field analysis predicts the threshold of the mixed phase (broken phase/ saturated phase) at 
/i ~ fi a . Future work should address possible relations between these observations. 

We believe that fi c indicates a physical critical point. All approaches (except the patho- 
logical quenched case) predict a transition or, at least, a clear change of behavior of ob- 
servables here. From the point of view of this section it is relatively easy to understand 
the robustness of this result : the probabilities plotted here underlie all the observables 
discussed earlier and the "flattness" of the distributions, which locates the critical point, is 
a qualitative feature which should appear in all the numerical procedures. 
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IV. OUTLOOK 



We have reached a partial understanding of the algorithms and we can point out some 
successes and failures. 

On the positive side, the method gives clear signs of the critical point \i c which should 
be the point of chiral symmetry restoration. The method also gives a current onset fi 
far different from \i c . This is most likely unphysical: it is not seen in the monomer-dimer 
results, it is the same as the pathological quenched onset, and it is the threshold of a phase 
characterized by forward propagation. Unfortunately, since Jo, and other observables such 
as the energy density, deviate from from zero at the unphysical fi Q point, their values near 
/i c cannot be trusted. So, although the present algorithm gives /i c accurately, it does not 
make any other phenomenologically reliable predictions. 

We expect that the early onset, fi Q , should disappear in a correct calculation. Physical 
arguments support this view as well as the monomer-dimer and strong coupling expansions 
discussed here. It might be that a high statistics run of the present algorithm will cancel 

In this case the method would be impractical, but, at least, not conceptually wrong. If 
this were true, we should develop a strategy to monitor the convergence of the method to 
the correct statistical ensemble, and to remove unphysical contributions to observables due 
to partial "cancellations" of unwanted onsets. 

A very unpleasant possibility, which we cannot exlude a priori, is that the results we are 
observing are indeed the final results at finite chemical potential with the Glasgow method. 
In this case, monomer-dimer simulations and strong coupling calculations would differ from 
the Glasgow results. This result would indicate intrinsic difficulties of the finite density 
lattice gauge theory simulation strategies. Some of these have been discussed in the text. 

We might also suspect that the problems stem from having generated configurations 
at zero chemical potential. This explanation is suggested from the standard problems en- 
countered by reweighting procedures, and from the behaviour observed in the Gross Neveu 
model || . In this case the Glasgow method can be improved if a better starting point were 
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invented. This is a worthwhile direction to pursue. 

At the present point, we have to accept that ensemble averaging does not help to suppress 
the pathologies of isolated configurations. It might well be that a satisfactory simulation of 
finite density QCD requires an algorithm which produces physical results on each configu- 
ration. A promising development of this sort is xQCD [ftm , where an irrelevant four fermi 
term is added to the standard QCD action used here. xQCD has the advantage that chiral 
symmetry breaking and the generation of a dynamical quark mass occurs configuration-by- 
configuration and the pion and sigma excitations are explicitly free of /i dependence. In fact, 
xQCD simulations do not suffer from the severe /i pathologies seen here [p6fl , but additional 
work, both theoretical and practical, is needed to see if xQCD really produces only physical 
results. Research in this topic is in progress. 
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FIGURES 
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FIG. 1. Quark number density from the Glasgow algorithm. m q = .1 on a 6 4 lattice. The onset 



/j, , and the saturation point (i s , are the same as the ones observed in the quenched approximation. 



The critical point for chiral symmetry restoration measured in a monomer-dimer calculation is 



/i c = .69(1), coincident with the little gap observed in our results. The same monomer-dimer 



results would, however, predict a very sharp transition with a critical density close to zero, in 



agreement with the results of the strong coupling expansion. 
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FIG. 3. Histogram of zeros accompanying Fig. |l|. Note a) the peak at [i c = .687(15) matching 
the small jump, to be contrasted with the monomer-dimer results fi c = .69(1). b) the corre- 
spondence of the extrema of the histogram with the onset [i of the current and its saturation 

Ms- 
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FIG. 4. Histogram of the zeros of the full model (top), and Histogram of the zeros of the 
determinant (bottom), hence of the quenched approximation. m q = .1, on a 8 4 lattice. 
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Re(exp(mu)) 

FIG. 5. Zeros in the e M plane for m = .1 8 4 lattice. The critical line is the thin line inside the 
denser region e M = e Mc . 
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FIG. 6. As in Fig. [|, but m q = .05, on a 6 4 lattice. 
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FIG. 7. As in Fig. |, but m q = .08, on a 6 4 lattice. 
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FIG. 9. Summary of the results for the critical point fi c , and current onset \x Q . \x c follows the 
prediction of the mean field analysis of ref.[18] (solid line). The onset is close to half the pion mass 
at small mass, and below half the pion mass for m q > .5 . 
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FIG. 10. Probability distributions for small chemical potentil at m 9 = .1 on the 8 4 lattice. The 
solid line is \i = 0, the dashed lines , from top to bottom at n = 0, are for \i = .1, .2, .3. 
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FIG. 11. Probability distributions around the onset fi for m q = .08 (top) and m q = .1 (bottom) 
on a 8 4 lattice. The leftmost histogram (solid) at m q = .08 is for \i = .28 , the rightmost is for 
H = .34. Bezier interpolations (from Gnuplot) are shown for fi = .28, .30, .32, .34 . At m q = .1 
fj, = .32, .34, .36, .38 from left to right. For both masses at fi Q the probability distribution moves on 
the positive n axes. 
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FIG. 12. Probability distributions in the critical region at m q = .1, on the 8 4 (top) and the 6 4 
(bottom), ijl is (.6, .683, .75) , from left to right (top) , and (.5. ,.6, .695, .75 , 1.) bottom. 
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FIG. 13. Probability distributions on the 8 lattice 



.1, for p = (.5 , .6 , .68, .683, .7 , .8 



, .9). Only the Bezier interpolations are shown. The complete results for several fj, values can be 



seen in Fig. 12 
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FIG. 14. Wo (solid) , and integrated probabilities W + and W (dash) at m 9 = .1, on the 8 4 
lattice 
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